On molecular and continuum boundary conditions for a miscible binary fluid 
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We show that molecular dynamics simulations can furnish useful boundary conditions at a solid 
surface bounding a two-component fluid. In contrast to some previous reports, convective-diffusive 
flow is consistent with continuum equations down to atomic scales. However, concentration gradients 
can produce flow without viscous dissipation that is inconsistent with the commonly used Navier 
slip condition. Also, differential wetting of the two components coupled to a concentration gradient 
can drive convective flows that could be used to make nano-pumps or motors. 
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Hydrodynamic theories of flow past a solid surface need 
to assume boundary conditions for the fluid velocities at 
the surface. Such interfacial behavior is often very diffi- 
cult to access experimentally. Recent simulation studies 
of fluids have revealed a range of boundary conditions for 
single component fluids and related them to the micro- 
scopic interactions fijfl]. However, there are many prob- 
lems where the appropriate boundary conditions are still 
in doubt. These include flow near a moving contact line, 
the liquid crystal order parameter in the presence of flow, 
and convective-diffusive flow of miscible fluids. 

In a recent letter |3| Koplik and Banavar studied two 
fluids that individually obeyed a "no-slip" boundary con- 
dition at the surface. In a mixture of these fluids, they 
found that the velocities of each component still vanished 
at the wall. They concluded that the boundary condi- 
tion accompanying the convective-diffusive transport of 
a binary fluid mixture demands that the velocities of the 
two species be equal at the wall. It was subsequently 
pointed out |J] that this boundary condition contradicts 
Fick's law of molecular diffusion, at least in its usual form 
where the diffusion coefficient is assumed to be position 
independent. Two possibilities were suggested. The first 
j| , that the convective flows studied in Ref . ||] were large 
enough to mask the effect of diffusion. The second, that 
the continuum equations need to be modified near the 
surface in order to obtain consistent results 

In an attempt to resolve this issue, Brenner and Gane- 
san H undertook a singular perturbation analysis. A 
continuum inner region with a refined form of Fick's law 
on which the no-slip condition applies for both compo- 
nents was asymptotically matched to an outer region 
which follows the standard form of Fick's law. As the 
details of the refined form of Fick's law are still some- 
what ambiguous, Brenner and Ganesan conclude that 
"simulations purporting to derive" boundary conditions 
"by direct probing of those molecules near to the wall 
will necessarily give rise to erroneous macroscale conclu- 
sions". Essentially, they argue that the simulations are 
useless since they do not give the information necessary 
to match up to the known equations for bulk fluids. If 



true, this would be a very disappointing situation. 

In this paper we examine flow boundary conditions for 
binary fluid mixtures using molecular dynamics simula- 
tions in the regime where diffusive flow is either domi- 
nant or on the same scale as the convective flow. We find 
bulk flow consistent with Fick's law within one molecular 
diameter from the wall, and appropriate boundary con- 
ditions for the known equations of bulk fluids are readily 
obtained. If concentration gradients are present, the in- 
dividual velocities of the two species are not equal at 
the wall. More surprisingly, we find that density gradi- 
ents can lead to a net flow without viscous dissipation! 
We also find that the coupling of wetting potentials with 
concentration gradients can give rise to Marangoni-like 
convective flow that could be used to make nano-pumps 
or nano-motors. 

We consider a mixture of two types of molecules, la- 
beled a and b, in a slit geometry similar to that used for 
a single fluid in Ref. jl). There are two walls, at z = 
and z = L, and periodic boundary conditions in the x — y 
plane. The walls contain type w atoms fixed to lattice 
sites forming two (111) layers of an fee surface. The inter- 
actions between atoms of type i and j separated by a dis- 
tance r are modeled using a Lennard-Jones (LJ) poten- 
tial, Vij(r) = idj [((Tij/r) 12 - (o-ij/r) 6 ] , where e# speci- 
fies the interaction energy and <jy the interaction length. 
Their averages are denoted by e and a, respectively. A 
characteristic time scale is given by r = o^m/e) 1 ' 2 , where 
m is the average of the molecular masses ruj. Unless 
specified, the force is truncated at r c — 2.2a to improve 
computational efficiency. 

To obtain a steady-state concentration gradient we ar- 
tificially construct a direction dependent osmotic mem- 
brane at x = 0. The membrane is less than one a thick, 
and preferentially transmits atoms of type a from left to 
right and conversely for atoms of type b. In a typical 
simulation we wait 5000r for the steady-state concentra- 
tion gradient to establish itself, and then average over a 
subsequent 15000 to 150000r to collect data. Note that 
these time scales are more than an order of magnitude 
greater than in Ref. jl, thus allowing much greater sen- 
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sitivity. Most channels studied had L — 16.4er, but we 
have also found similar results for systems twice as wide. 
The diffusive flux Jj of particles of type j, 

J i = Pj( v j ~ v), (1) 

is defined Q relative to the barycentric, or mass aver- 
aged, velocity v = piVi)/p, where pi is the position- 
dependent mass density of species i, and p = - pj is the 
total mass density. Fick's law relates the diffusive flux to 
the gradient of the concentration Cj = Pj/p (?]] , 

3 3 =-pDVc 3 , (2) 

where D is the diffusion constant. The results of Refs. 

suggest that there may be a macroscopic layer where 
the diffusion constant depends on proximity to the wall. 

For simplicity, we first consider a system where ov, = a 
for all interactions and m a — mi, — m. All fluid-fluid 
€ij = e, and the wall-fluid coupling e w f = e wa — e w b is 
varied. In this case entropy is the only driving force. Fig- 
ure |](a) shows the resulting steady state densities pj of 
the two species. Both vary linearly and the total density 
is constant throughout the system. Since a and b parti- 
cles have equal and opposite concentration gradients we 
find that their fluxes are equal and opposite and v = 0. 
If the difference between individual velocities vanished at 
the wall, as in Ref. 4, both fluxes would have to vanish 
there. Figure [j](b) illustrates that the fluxes may actually 
increase at the wall. To eliminate variations with x we 
plot D = v x ,j(cj /d x Cj) as a function of height z for three 
different wall-fluid interaction strengths. The diffusion 
constant, and thus the flux, goes down near the wall for 
£u;//e = 1-0, remains constant for e w f/e — 0.3, and goes 
up for e w f/e = 0.1. Thus, for e w f/e = 0.1, the difference 
between the individual velocities of the two species is ac- 
tually highest at the wall. The value of D is normalized 
by the bulk diffusion constant which was calculated for a 
system where the walls were replaced by periodic bound- 
ary conditions in the z direction. As can be seen, there is 
no statistical difference from bulk diffusion at distances 
more than about one molecular diameter from the wall. 
Thus, in contrast to References f|-^|, there is no need to 
supplement Fick's law with a continuum boundary layer 
which obeys a modified constitutive relation. 

Earlier work jl],|j on single component LJ fluids is well 
described by the Navier boundary condition, which as- 
sumes that the velocity difference between solid and fluid 
is proportional to the viscous shear stress. For our ge- 
ometry this implies 

Vx\wall = L s d z V x \ wa ll, (3) 

where L s is called the slip length. A larger value of L s 
implies less drag at the interface, and we find that it cor- 
relates with greater diffusion at the wall. In most cases, 
the slip length has atomic dimensions (l] Q , and it is less 
than 2a for single fluids with the walls considered here. 



The Navier condition presupposes that in the purely 
diffusive case, where d z v x — 0, the barycentric velocity 
vanishes identically. However, Fig. ^(a) shows that this is 
not necessarily the case. For this simulation the masses 
were changed to m a = 0.75m and nib — 1.25m and 
e wf /t — 1-0. We find that the number densities Ni are the 
same, within statistical errors, as those of Fig. [l](a) and 
the mean velocities of each species are also unchanged. 
However, due to the mass difference, there is now a mass 
density gradient and a small net mass flux J = pv along 
the channel (Fig. |^(a)). One might wonder if the Navier 
condition still applies, but with a very small velocity gra- 
dient and very large slip length. However, both single and 
two-component simulations of purely convective flow for 
the same parameters yield "stick" boundary conditions, 
i.e. the effective "wall" position is inside the fluid and 
L s is essentially zero. Moreover, there is no measurable 
viscous stress and the entire flux profile follows from a 
general continuum relation that we now describe. 

For the example of Fig. |^(a) the mean molar velocity 
v m = (J^i NiVi)/N vanishes, where Ni is the molar den- 
sity of component i and N = Y^a Aj- Using the definition 
of the diffusive flux (Eq. (§)), and Fick's law (Eq. (§)) 
one can show that 

v = v m + DV [HN/p)\ . (4) 

Thus the barycentric and molar velocities differ when 
there is a diffusive flow driven by a gradient in a to- 
tal density (N or p). In the case considered in Fig. [l], 
the total densities are constant and both v and v m van- 
ish everywhere for purely diffusive flow. For the case 
of Fig. ||(a), v m vanishes and N is essentially constant. 
Thus Eq. |] implies J = pv = —DV p. The dashed line in 
Fig. ||(a) shows that this expression reproduces the ob- 
served mass flux. Note that the height dependence comes 
entirely from the variation in D that is also evident in 
Fig. [j](b). If smaller values of e w f/e are used, the flux 
may actually increase at the wall. This leads to different 
signs on the two sides of Eq. || 

The observed failure of the Navier condition has im- 
portant consequences. Numerous arguments related to 
the theory of diffusive flow, such as those in Ref. B , rely 
on the Navier condition to enforce a vanishing barycen- 
tric velocity in the absence of diffusion. The conclusions 
of these works need to be reexamined. One may won- 
der whether the correct boundary condition for purely 
diffusive flow is the identical vanishing of some other 
mean velocity such as v m . However, it is easy to con- 
struct exceptions to such ansatzes. The difficulty with 
such boundary conditions is most readily seen by consid- 
ering steady state situations. Steady state implies that 
dtp = dfNi = and the corresponding continuity equa- 
tions then imply that 

d a (pv a ) = d a (Nv™)=0. (5) 
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Multiplying Eq. (Q) by N and taking the divergence 
yields 

d a (Nv a ) = d a (Nv™) + d a (NDd a [HN/p)}) . (6) 

Rewriting d a (Nv a ) as d a [(N/p)(pv a )} and making use 
of Eq. (||) , one can easily obtain 

v a d a (N/p) = -d a (NDd a {HN/p)}) . (7) 
P 

An analogous relation for the mean molar velocity can 
easily be derived, 

v™d a (p/N) = ±d a (pDd a [HN/p)}) . (8) 

Thus the continuum relations themselves place con- 
straints on mean velocities. Any additional requirement 
that a mean velocity vanish imposes constraints on the 
gradients of the density. While these may be consistent 
with some situations, they will in general over-specify the 
problem. 

We have tested the above equations with simulations 
using a range of atom sizes and masses. To illustrate 
their application, we use the parameters for Fig. ||(a), 
but adjust the probability of reflection at the osmotic 
wall to produce different net average mass flows down 
the channel. Typical examples of the resulting steady 
state mass density and flow profiles for e w f/c = 0.3 are 
shown in Fig. ^(b) and (c). Neither the barycentric nor 
the mean molar velocity vanishes. Indeed, they are in- 
dependent of height at distances more than a molecu- 
lar diameter from the wall. We find that the mean mo- 
lar density N is independent of x, and that only the x 
components of the velocity are nonzero. Thus Eq. (|^) 
simplifies to v" 1 — Dd x p/d x p. The spatially averaged 
versus d x p/d x p for several different osmotic bound- 
ary conditions is shown in Fig. ||(d). We find a good fit 
to a straight line whose slope D — 0.067 ± 0.004 agrees 
with the value D = 0.070 ± 0.004 found directly from 
Fick's law (Eq. (|)). 

It is clear from the above examples that for multi- 
component fluids the Navier condition must be modi- 
fied by subtracting the diffusive flow from the left hand 
side of Eq.(||). We find that simulations of a variety of 
convective-diffusive flows can then be fit with a common 
value of the slip length || . 

So far we have examined only neutral wetting where 
both a and b have the same interaction with the walls. 
In many realistic cases, one component will preferentially 
wet the wall. To examine this situation we return to the 
case where the two particles are indistinguishable except 
for their labels. For particles a (b) we now change the 
wall-fluid interaction at the top (bottom) wall so that 
it is purely repulsive by truncating the potential at its 
minimum rather than at 2.2a. Recall that previously 
for this fluid we found that all average velocities were 



zero. Figure [|(a) shows the average velocity produced 
by changing the wetting properties. One sees a remark- 
ably strong shear flow (~m/s) with significant "slip" at 
the stationary walls. This slip has the wrong sign for the 
Navier condition and its magnitude is inconsistent with 
results for convective flows. 

The driving force for the flow in Fig. |^(a) comes from 
the externally imposed concentration gradient and the 
variation of the interfacial free energy of the walls with 
concentration. Figure ^|(b) shows the mass fraction c a 
at different cross-sections along the channel. As can be 
seen from the figure, along the walls the system prefers 
the more strongly wetting fluid species. The difference 
between the concentration at the wall and in the center 
of the channel increases with x. This increase leads to a 
rise in interfacial tension. Values of the surface tension 
7, calculated using the mechanical definition of Kirkwood 
and Buff JIo| , are shown in Fig. |^(c) . The boundary con- 
dition relating the shear stress on the wall a xz \ w to the 
viscous stress in the fluid <J xz \f is analogous to that used 
for Marangoni flow at a two fluid boundary []l]J : 

<Txz\w ~ &xz\f = d x j. (9) 

It can be shown [^) that there is also a velocity discon- 
tinuity given by an integral of the Marangoni stress over 
the interfacial region that produces the slip in Fig. ||(a). 
This can be incorporated into a generalized Navier slip 
condition by adding a source term proportional to d x ^. 
These boundary conditions determine the net stress on 
the wall. The result is consistent with stresses measured 
directly in our simulations, and ranges from 50 to 75% 
of <9 X 7 (^MPa) for the cases we have studied. This force 
could be used to drive a nano-motor. We have also con- 
structed systems where both the top and bottom walls 
prefer the same component. In this case, one generates 
Poiseuille-like flow which could be used in a nano-pump. 

In conclusion, we find that the flow boundary condi- 
tion for convective-diffusive flow is not one of equal ve- 
locities for all species. In addition, average velocities do 
not, in general, vanish at the wall in the absence of vis- 
cous stress. Diffusive mass transport can contribute to 
a significant average velocity at the wall. Further, the 
presence of concentration gradients along the wall in the 
general case of non-neutral wetting can result in signifi- 
cant Marangoni- type forces which drive convective flow. 
These effects should be readily applicable to the design 
of new micro-fluidic devices and may be relevant to the 
function of numerous biological systems. 

This material is based upon work supported by the 
National Science Foundation under Grant No. 0083286. 
We thank Intel Corporation for donating the worksta- 
tions used for our simulations, which were performed us- 
ing LAMMPS from Sandia National Laboratories. 
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FIG. 2. (a) The observed mass flux J = pv x (♦) and that 
predicted J = -Dd x p (*) by Eq. @ for v m = 0. (b) Mass 
density of fluid particles of type a (♦) and 6 (*) and total den- 
sity (■). (c) The mass flux normalized by the mean density 
po, v°z V = pVx/po- (d) Mean molar velocity versus d^p/d x p 
for different osmotic walls. The line is a linear fit giving 
D = 0.067 ±0.004cr 2 /V- 
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FIG. 1. (a) Density of fluid particles of type a (♦) and 
b (*), and total density (■). The osmotic membrane is lo- 
cated at x = 0. The periodic cell dimensions along the x 
and y-directions are 42.5cr and 12.3<r, respectively, (b) Diffu- 
sion constant for e w f/e = 1.0 (■), 0.3 (*), and 0.1 (♦). We 
normalize by the bulk diffusion constant, Db = 0.068<t 2 /t, 
measured in a system without walls. Away from the osmotic 
membrane, the value of D is independent of x and y. 
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FIG. 3. (a) Barycentric velocity down the channel for a 
case where the wall at z = preferentially wets a and the 
other wall preferentially wets b. (b) Mass fraction c a — p a /p 
versus height evaluated at planes 1/4 (a) and 3/4 (★) of the 
way along the channel in the a;-direction. (c) Surface tension 
along the top (♦) and bottom (*) walls. 
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